clear; close all; clc;

%% load results

cdir = pwd;
savef = 1;                  % savef = 1 to save plots
fdir  = 'Figures/Crisis';  % location to save plots

cname = 'Argentina';
odir  = strcat(cname,'/OUTPUT_1');  % benchmark calibration for Argentina

cd(odir)

%--grids
load bvec.txt;
load zvec.txt;
load Pzvec.txt;

%--policies
load bpol.txt;
load dpol.txt;
load bopt_loc.txt;
load Rsched.txt;

load prms.txt;

%--simulation
load bsim.txt;
load cnsim.txt;
load idsim.txt;
load rhosim.txt;
load nsim.txt;
load gsim.txt;
load zsim.txt;
load issim.txt;
load ysim.txt;
load qbsim.txt;
load qsim.txt;
load psim.txt;

cd(cdir)

%prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
Rs  = prms(1) ; dlta = prms(2) ; ybarL = prms(3) ; ybarH = prms(4) ; kppa  = prms(5);
pgL = prms(6) ; pgH  = prms(7) ; bta   = prms(8) ; gL    = prms(9) ; gH    = prms(10); 
sdz = prms(11); psB   = prms(12);

psG   = 1.0-psB;

gvec = [gL, gH];
Pg   = [pgL  , 1-pgL; ...
        1-pgH, pgH];
Ps   = [psB  , 1-psB; ...
        1-psG, psG];

Nb = numel(bvec); Nz = numel(zvec); Ng = 2; Ns = 2;


%% simulation computations

Tsim  = numel(ysim);
Tburn = 1000;
TTvec = (1:1:Tsim)';

%---adjust variables
idsim = idsim(2:Tsim+1);

rhosimx = rhosim;
rhosim  = rhosimx - (Rs-1);

nsim  = nsim./ysim;
bsim  = bsim(1:Tsim)./ysim; 
fsim  = cnsim(1:Tsim)./ysim; 
qbsim = qbsim(1:Tsim)./ysim; qbsim = qbsim/dlta;

dnsim  = nsim(2:Tsim) -nsim(1:Tsim-1) ; dnsim  = [0; dnsim];
dbsim  = bsim(2:Tsim) -bsim(1:Tsim-1) ; dbsim  = [0; dbsim];
dfsim  = fsim(2:Tsim) -fsim(1:Tsim-1) ; dfsim  = [0; dfsim];
dqbsim = qbsim(2:Tsim)-qbsim(1:Tsim-1); dwbsim = [0; dqbsim];
dysim = log(ysim(2:Tsim))-log(ysim(1:Tsim-1));  dysim = [0; dysim];

%---set condition
rhosimax = rhosim(Tburn+1:Tsim); idsimx = idsim(Tburn+1:Tsim);

rhoxnd  = rhosimax(idsimx == 0);

rhoxmn  = mean(rhoxnd);
rhoxsd  = std(rhoxnd);
rhocond = rhoxmn + 1.0*rhoxsd;

xcond = 0*ones(Tsim,1);
for tt = 2+1:Tsim-1-1    
    %if (xcond(tt-1) == 0 && idsim(tt-1) == 0 && idsim(tt) == 0 && rhosim(tt) > rhosim(tt-1) && rhosim(tt) > rhosim( tt+1) && rhosim(tt) > rhosim(tt+2) && rhosim(tt) > rhocond)
    if (xcond(tt-1) == 0 && idsim(tt-1) == 0 && idsim(tt) == 0 && rhosim(tt) > rhosim(tt-2) && rhosim(tt) > rhosim(tt-1) && rhosim(tt) > rhosim(tt+1) && rhosim(tt) > rhosim(tt+2) && rhosim(tt) > rhocond)    
        xcond(tt) = 1;        
    end
end

xcond = xcond.*(TTvec > Tburn);
tcond = TTvec(xcond == 1);

Ncond = sum(xcond);
LAG  = 3;
FRW  = 3;
NIRF = LAG + 1 + FRW;
IRF_xc_r  = 189*ones(Ncond,NIRF);  IRF_r  = 189*ones(NIRF,1);  IRF_sd_r  = 189*ones(NIRF,1);  IRF_xs_r  = 189*ones(Ns,NIRF); IRF_xs_sd_r  = 189*ones(Ns,NIRF);
IRF_xc_n  = 189*ones(Ncond,NIRF);  IRF_n  = 189*ones(NIRF,1);  IRF_sd_n  = 189*ones(NIRF,1);  IRF_xs_n  = 189*ones(Ns,NIRF); IRF_xs_sd_n  = 189*ones(Ns,NIRF);
IRF_xc_b  = 189*ones(Ncond,NIRF);  IRF_b  = 189*ones(NIRF,1);  IRF_sd_b  = 189*ones(NIRF,1);  IRF_xs_b  = 189*ones(Ns,NIRF); IRF_xs_sd_b  = 189*ones(Ns,NIRF);
IRF_xc_f  = 189*ones(Ncond,NIRF);  IRF_f  = 189*ones(NIRF,1);  IRF_sd_f  = 189*ones(NIRF,1);  IRF_xs_f  = 189*ones(Ns,NIRF); IRF_xs_sd_f  = 189*ones(Ns,NIRF);
IRF_xc_qb = 189*ones(Ncond,NIRF);  IRF_qb = 189*ones(NIRF,1);  IRF_sd_qb = 189*ones(NIRF,1);  IRF_xs_qb = 189*ones(Ns,NIRF); IRF_xs_sd_qb = 189*ones(Ns,NIRF);
IRF_xc_dn = 189*ones(Ncond,NIRF);  IRF_dn = 189*ones(NIRF,1);  IRF_sd_dn = 189*ones(NIRF,1);  IRF_xs_dn = 189*ones(Ns,NIRF); IRF_xs_sd_dn = 189*ones(Ns,NIRF);
IRF_xc_db = 189*ones(Ncond,NIRF);  IRF_db = 189*ones(NIRF,1);  IRF_sd_db = 189*ones(NIRF,1);  IRF_xs_db = 189*ones(Ns,NIRF); IRF_xs_sd_db = 189*ones(Ns,NIRF);
IRF_xc_df = 189*ones(Ncond,NIRF);  IRF_df = 189*ones(NIRF,1);  IRF_sd_df = 189*ones(NIRF,1);  IRF_xs_df = 189*ones(Ns,NIRF); IRF_xs_sd_df = 189*ones(Ns,NIRF);
IRF_xc_s  = 189*ones(Ncond,NIRF);  IRF_s  = 189*ones(NIRF,1);  %IRF_sd_b  = 189*ones(NIRF,1);  
IRF_xc_y  = 189*ones(Ncond,NIRF);  IRF_y  = 189*ones(NIRF,1);  IRF_sd_y  = 189*ones(NIRF,1);  
IRF_xc_dy = 189*ones(Ncond,NIRF);  IRF_dy = 189*ones(NIRF,1);  IRF_sd_dy = 189*ones(NIRF,1);  

IRF_p05_r  = 189*ones(NIRF,1); IRF_p95_r  = 189*ones(NIRF,1);
IRF_p05_qb = 189*ones(NIRF,1); IRF_p95_qb = 189*ones(NIRF,1);

dw = 4;
for icond = 1:Ncond
    disp(['icond = ',num2str(icond),'/',num2str(Ncond)])
    tt    = tcond(icond);
    ttind = (tt-LAG:tt+FRW)';
    
    IRF_xc_r(icond,:)  = rhosim(ttind);
    IRF_xc_n(icond,:)  = nsim(ttind);
    IRF_xc_b(icond,:)  = bsim(ttind);
    IRF_xc_qb(icond,:) = qbsim(ttind);
    IRF_xc_f(icond,:)  = fsim(ttind);
    IRF_xc_dn(icond,:) = dnsim(ttind);
    IRF_xc_db(icond,:) = dbsim(ttind); 
    IRF_xc_df(icond,:) = dfsim(ttind); 
    IRF_xc_y(icond,:)  = ysim(ttind); 
    IRF_xc_dy(icond,:) = dysim(ttind); 

    IRF_xc_s(icond,:) = issim(ttind); 
    
end

llx = LAG+1; sdata  = IRF_xc_s(:,llx); 

i10 = floor(0.05*Ncond); i90 = ceil(0.95*Ncond);
for ll = 1:NIRF
    IRF_r(ll)  = mean(IRF_xc_r(:,ll));  IRF_sd_r(ll)  = std(IRF_xc_r(:,ll));
    IRF_n(ll)  = mean(IRF_xc_n(:,ll));  IRF_sd_n(ll)  = std(IRF_xc_n(:,ll));
    IRF_b(ll)  = mean(IRF_xc_b(:,ll));  IRF_sd_b(ll)  = std(IRF_xc_b(:,ll));
    IRF_f(ll)  = mean(IRF_xc_f(:,ll));  IRF_sd_f(ll)  = std(IRF_xc_f(:,ll));
    IRF_qb(ll) = mean(IRF_xc_qb(:,ll)); IRF_sd_qb(ll) = std(IRF_xc_qb(:,ll));
    IRF_dn(ll) = mean(IRF_xc_dn(:,ll)); IRF_sd_dn(ll) = std(IRF_xc_dn(:,ll)); 
    IRF_db(ll) = mean(IRF_xc_db(:,ll)); IRF_sd_db(ll) = std(IRF_xc_db(:,ll)); 
    IRF_df(ll) = mean(IRF_xc_df(:,ll)); IRF_sd_df(ll) = std(IRF_xc_df(:,ll)); 

    yy = IRF_xc_r(:,ll) ; yys = sort(yy); IRF_p05_r(ll)  = yys(i10); IRF_p95_r(ll)  = yys(i90); 
    yy = IRF_xc_qb(:,ll); yys = sort(yy); IRF_p05_qb(ll) = yys(i10); IRF_p95_qb(ll) = yys(i90); 

    for is = 1:Ns       
       ydatax = IRF_xc_r(:,ll) ; ydata = ydatax(sdata==is); IRF_xs_r(is,ll)  = mean(ydata); IRF_xs_sd_r(is,ll)  = std(ydata);
       ydatax = IRF_xc_n(:,ll) ; ydata = ydatax(sdata==is); IRF_xs_n(is,ll)  = mean(ydata); IRF_xs_sd_n(is,ll)  = std(ydata);
       ydatax = IRF_xc_b(:,ll) ; ydata = ydatax(sdata==is); IRF_xs_b(is,ll)  = mean(ydata); IRF_xs_sd_b(is,ll)  = std(ydata);
       ydatax = IRF_xc_f(:,ll) ; ydata = ydatax(sdata==is); IRF_xs_f(is,ll)  = mean(ydata); IRF_xs_sd_f(is,ll)  = std(ydata);

       ydatax = IRF_xc_dn(:,ll); ydata = ydatax(sdata==is); IRF_xs_dn(is,ll) = mean(ydata); IRF_xs_sd_dn(is,ll) = std(ydata);
       ydatax = IRF_xc_db(:,ll); ydata = ydatax(sdata==is); IRF_xs_db(is,ll) = mean(ydata); IRF_xs_sd_db(is,ll) = std(ydata);
       ydatax = IRF_xc_df(:,ll); ydata = ydatax(sdata==is); IRF_xs_df(is,ll) = mean(ydata); IRF_xs_sd_df(is,ll) = std(ydata);

    end
end

IRF_lb_r  = IRF_r -IRF_sd_r ; IRF_ub_r  = IRF_r +IRF_sd_r;  IRF_xs_lb_r  = IRF_xs_r -IRF_xs_sd_r ; IRF_xs_ub_r  = IRF_xs_r +IRF_xs_sd_r ;
IRF_lb_n  = IRF_n -IRF_sd_n ; IRF_ub_n  = IRF_n +IRF_sd_n;  IRF_xs_lb_n  = IRF_xs_n -IRF_xs_sd_n ; IRF_xs_ub_n  = IRF_xs_n +IRF_xs_sd_n ;
IRF_lb_b  = IRF_b -IRF_sd_b ; IRF_ub_b  = IRF_b +IRF_sd_b;  IRF_xs_lb_b  = IRF_xs_b -IRF_xs_sd_b ; IRF_xs_ub_b  = IRF_xs_b +IRF_xs_sd_b ;
IRF_lb_f  = IRF_f -IRF_sd_f ; IRF_ub_f  = IRF_f +IRF_sd_f;  IRF_xs_lb_f  = IRF_xs_f -IRF_xs_sd_f ; IRF_xs_ub_f  = IRF_xs_f +IRF_xs_sd_f ;
IRF_lb_qb = IRF_qb-IRF_sd_qb; IRF_ub_qb = IRF_qb+IRF_sd_qb; 
IRF_lb_dn = IRF_dn-IRF_sd_dn; IRF_ub_dn = IRF_dn+IRF_sd_dn; IRF_xs_lb_dn = IRF_xs_dn-IRF_xs_sd_dn; IRF_xs_ub_dn = IRF_xs_dn+IRF_xs_sd_dn;
IRF_lb_db = IRF_db-IRF_sd_db; IRF_ub_db = IRF_db+IRF_sd_db; IRF_xs_lb_db = IRF_xs_db-IRF_xs_sd_db; IRF_xs_ub_db = IRF_xs_db+IRF_xs_sd_db;
IRF_lb_df = IRF_df-IRF_sd_df; IRF_ub_df = IRF_df+IRF_sd_df; IRF_xs_lb_df = IRF_xs_df-IRF_xs_sd_df; IRF_xs_ub_df = IRF_xs_df+IRF_xs_sd_df;



%% IRF + confidence band plot

linesz   = 0.5;
labelsz  = 9;
legendsz = 9;
numsz    = 9;
textsz   = 9;

IRF_xaxis = -LAG:1:FRW; IRF_xaxis = IRF_xaxis'; IRF_T = numel(IRF_xaxis);
xx = [IRF_xaxis(1:IRF_T)', IRF_xaxis(IRF_T:-1:1)'];
rrlb = IRF_lb_r' ; rrub = IRF_ub_r' ; YYrr = [rrlb, rrub(IRF_T:-1:1)];
bblb = IRF_lb_b' ; bbub = IRF_ub_b' ; YYbb = [bblb, bbub(IRF_T:-1:1)];
fflb = IRF_lb_f' ; ffub = IRF_ub_f' ; YYff = [fflb, ffub(IRF_T:-1:1)];
nnlb = IRF_lb_n' ; nnub = IRF_ub_n' ; YYnn = [nnlb, nnub(IRF_T:-1:1)];
qblb = IRF_lb_qb'; qbub = IRF_ub_qb'; YYqb = [qblb, qbub(IRF_T:-1:1)];

rrlb = IRF_p05_r' ; rrub = IRF_p95_r' ; YYpprr = [rrlb, rrub(IRF_T:-1:1)];
qblb = IRF_p05_qb'; qbub = IRF_p95_qb'; YYppqb = [qblb, qbub(IRF_T:-1:1)];

CY   = [0.8 0.3 0.3];
Cirf = [0.4 0.4 0.4];


fig = figure(5000+1); clf
fill(xx,100*YYpprr,CY,'edgecolor',CY,'FaceAlpha',0.50); hold on
plot(IRF_xaxis,100*IRF_r,'color',Cirf,'LineWidth',linesz)
set(gca,'XGrid','off','YGrid','off','Fontsize',numsz)
xlabel('Periods afer debt crisis','Interpreter','LaTex','Fontsize',labelsz)
ylabel('spread $\rho - r^*$ (\%)','Interpreter','LaTex','Fontsize',labelsz)
hold off
if savef == 1    
    paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
    paper_size = [3.0 2.0]; % [width height]
    paper_position = [0 0 3.0 2.0]; % [left bottom width height].
    paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
    fig = gcf;
    fig.PaperUnits = paper_unit;
    fig.PaperSize = paper_size;
    fig.PaperPosition = paper_position;
    fig.PaperOrientation = paper_orientation;
    cd(fdir)
    saveas(gcf,strcat('crisis_IRF_r_',cname,'.pdf'))    
    cd(cdir)    
end

fig = figure(5000+2); clf
fill(xx,100*YYppqb,CY,'edgecolor',CY,'FaceAlpha',0.50); hold on
plot(IRF_xaxis,100*IRF_qb,'color',Cirf,'LineWidth',linesz)
set(gca,'XGrid','off','YGrid','off','Fontsize',numsz)
xlabel('Periods afer debt crisis','Interpreter','LaTex','Fontsize',labelsz)
ylabel('debt market value(\%)','Interpreter','LaTex','Fontsize',labelsz)
hold off
if savef == 1
    %run graph_extended.m
    paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
    paper_size = [3.0 2.0]; % [width height]
    paper_position = [0 0 3.0 2.0]; % [left bottom width height].
    paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
    fig = gcf;
    fig.PaperUnits = paper_unit;
    fig.PaperSize = paper_size;
    fig.PaperPosition = paper_position;
    fig.PaperOrientation = paper_orientation;
    cd(fdir)
    saveas(gcf,strcat('crisis_IRF_qb_',cname,'.pdf'))    
    cd(cdir)    
end


